function [Yp,numsv] = mmkp(V,Y,C)
tic
l = size(Y,1);
[n,s] = size(V);

cls = unique(Y);
c = size(cls,1);
nump = c*(c-1)/2;

sz = zeros(c,1);
idx = cell(c,1);
for i = 1:c
    idx{i}=find(Y==cls(i));
    sz(i)=size(idx{i},1);
end

ix = cell(c,c);
Y = cell(c,c);
P = cell(c,c);

for i = 1:c
    for j = i+1:c
        ix{i,j} = [idx{i};idx{j}];
        Y{i,j} = [ones(sz(i),1);-ones(sz(j),1)];
        P{i,j} = repmat(Y{i,j},1,s).*V(ix{i,j},:);
    end
end

disp('Start optimization...');
tic
ct = 1;
if(nargin == 2)
    cvx_begin sdp quiet
    cvx_solver sedumi
    cvx_solver_settings('maxiter',500);
    variable t
    variable gam(s,nump)
    variable M(s,s) symmetric
    variable bb(nump,1)
    minimize (t);
    for i = 1:c
        for j = i+1:c
            [M gam(:,ct);gam(:,ct)' t] >= 0;
            P{i,j}*gam(:,ct)+bb(ct)*Y{i,j}-ones(size(P{i,j},1),1) >= 0;
            ct=ct+1;
        end
    end
    trace(M) <= s;
    cvx_end
else
    cvx_begin sdp quiet
    cvx_solver sedumi
    cvx_solver_settings('maxiter',500);
    variable t
    variable gam(s,nump)
    variable M(s,s) symmetric
    variable del(l,1) nonnegative
    variable bb(nump,1)
    minimize (t+2*C*sum(del));
    for i = 1:c
        for j = i+1:c
            [M gam(:,ct);gam(:,ct)' t] >= 0;
            P{i,j}*gam(:,ct)+bb(ct)*Y{i,j}-ones(size(P{i,j},1),1)+del(ix{i,j}) >= 0;
            ct=ct+1;
        end
    end
    trace(M) <= s;
    cvx_end
end
optime = toc

VM = V*gam+repmat(bb',n,1);
vot = zeros(n,c);

ct = 1;
for i = 1:c
    for j = i+1:c
        vot(:,i) = vot(:,i)+(VM(:,ct)>0);
        vot(:,j) = vot(:,j)+(VM(:,ct)<0);
        ct=ct+1;
    end
end
[~,Yp] = max(vot,[],2);
if (nargin == 2)
    numsv = 0;
else
    numsv = sum(del>0);
end
Yp = Yp-1;
end